epi draft 



Effects of an embedding bulk fluid on phase separation dynamics 
in a thin liquid film 



o 



O 



S. Ramachandran^ , S. KoMURA^ and G. Gompper^ 

^ Department of Chemistry, Graduate School of Science and Engineering, Tokyo Metropolitan University, Tokyo 192- 
)0397, Japan. 

Institut fiir Festkorperforsehung, Forschungszentrum Jiilich - D-52425 Jiilich, Germany, EU 



PACS 68 . 05 . -n - Liquid-liquid interfaces 

PACS 87. 16.D — Membranes, bilayers, and vesicles 

PACS 87. 16. dp - Transport, including channels, pores. 



and lateral difFusion 



Abstract. - Using dissipative particle dynamics simulations, we study the efltects of an embedding 
bulk fluid on the phase separation dynamics in a thin planar liquid film. The domain growth 
exponent is altered from 2D to 3D behavior upon the addition of a bulk fluid, even though the 
phase separation occurs in 2D geometry. Correlated diffusion measurements in the film show that 
the presence of bulk fluid changes the nature of the longitudinal coupling diffusion coefficient from 
logarithmic to algebraic dependence of 1/s, where s is the distance between the two particles. This 
result, along with the scaling exponents, suggests that the phase separation takes place through 
the Brownian coagulation process. 
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, Introduction. — Lipid molecules constituting the 
ihembranes of biological cells play a major role in the reg- 
ulation of various cellular processes. About a decade ago 
Simons and Ikonen proposed a hypothesis which suggests 
that the lipids organize themselves into sub-micron sized 
domains termed "rafts" [T]. The rafts serve as platforms 
for proteins, which in turn attributes a certain function- 
ality to each domain. Although there have been extensive 
studies in this area, the details of the underlying physical 
itiechanisms leading to formation of rafts, their stability, 
and the regulation of their finite domain size remain elu- 
sive. Numerous experiments on intact cells and artificial 
membranes containing mixtures of saturated lipids, un- 
saturated lipids and cholesterol, have demonstrated the 
segregation of the lipids into liquid-ordered and liquid- 
disordered phases [5]. Recent experimental observations 
of critical fluctuations point towards the idea that the 
cell maintains its membrane slightly above the critical 
point [3113]. Below the transition temperature, there have 
been studies on the dynamics in multicomponent mem- 
branes such as diffusion of domains [5] and domain coars- 
ening [S]. A clear understanding of phase separation may 
contribute towards a better explanation of the dynamics 
of lipid organization in cell membranes. Apart from bi- 
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ological membranes, it is of relevance to understand the 
dynamics of Langmuir monolayer systems which are also 
thin fluid films. 

Phase separation of binary fluids following a quench 
has been under study for over forty years 0- The dy- 
namic scaling hypothesis assumes that there exists a scal- 
ing regime characterized by the average domain size R 
that grows with time t as R ^ t"' with an universal ex- 
ponent a. For three-dimensional (3D) off-critical binary 
fluids, there is an initial growth by the Brownian coagu- 
lation process [8], followed by the Lifshitz-Slyozov (LS) 
evaporation-condensation process [3]; both mechanisms 
show a growth exponent a = 1/3. This is followed by 
a late time inertial regime of a = 2/3 [TO]. For critical 
mixtures, there is an intermediate a = I regime owing 
to interface diffusion [11]. The scenario is slightly differ- 
ent for pure two-dimensional (2D) systems |12| . For an 
off-critical mixture, it was predicted that after the initial 
formation of domains, they grow by the Brownian coagu- 
lation mechanism with a different exponent a = 1/2 (as 
will be explained later), followed by a crossover to the 
LS mechanism which gives a = 1/3 even in 2D. For crit- 
ical mixtures, on the other hand, the initial quench pro- 
duces an interconnected structure which coarsens and then 
breaks up due to the interface diffusion with an exponent 
a = 1/2. After the breakup processes, coarsening takes 
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place through Brownian coagulation that is again charac- 
terized by the a = 1/2 scaling [5]. These predictions were 
confirmed by molecular dynamics simulations in 2D |13j . 
The exponent a = 1/2 was also observed in 2D lattice- 
Boltzmann simulations in the presence of thermal noise 
for a critical mixture [14] . 

Although biomembranes composed of lipid bilayers can 
be regarded as 2D viscous fluids, they are not isolated 
pure 2D systems since lipids are coupled to the adjacent 
fluid. Hence it is of great interest to investigate the phase 
separation dynamics in such a quasi-2D liquid film in the 
presence of hydrodynamic interaction. (We use the word 
"quasi-2D" whenever the film is coupled to the bulk fluid.) 
To address this problem, wc consider a 2D binary viscous 
fluid in contact with a bulk fluid. Our approximation of 
the membrane as a planar 2D liquid film is valid in the 
limit of large bending rigidity (common in biological mem- 
branes) or in the presence of a lateral tension, which both 
act to suppress membrane undulations. We employ a sim- 
ple model in which the film is confined to a plane with 
the bulk fluid particles added above and below. In our 
model using dissipative particle dynamics (DPD) simula- 
tion technique, the exchange of momentum between the 
film and the bulk fluid is naturally taken into account. 
We particularly focus on the effect of bulk fluid on the 
quasi-2D phase separation. Wc show that the presence of a 
bulk fluid will alter the domain growth exponent from that 
of 2D to 3D indicating the significant role played by the 
film-solvent coupling. In order to elucidate the underlying 
physical mechanism of this effect, we have looked into the 
diffusion properties in the film by measuring two-particle 
correlated diffusion. Our result suggests that quasi-2D 
phase separation proceeds by the Brownian coagulation 
mechanism which reflects the 3D nature of the bulk fluid. 
Such a behavior is universal as long as the domain size 
exceeds the Saffman-Delbriick length [15] . 



Model and simulation technique. For the pur- 
pose of our study, we use a structureless model of the 2D 
liquid film within the DPD framework (TBIITJ. As shown 
in fig. [U the 2D film is represented by a single layer of par- 
ticles confined to a plane. In order to study phase separa- 
tion, we introduce two species of particles, A and B. The 
bulk fluid which wc call as "solvent" [S) is also represented 
by single particles of same size as that of the film particles. 
All particles have the same mass m. We avoid using the 
existing DPD models for a self-assembling bilayer P^[T5] 
as they inherently include bending and protrusion modes, 
which makes it difhcult to separate hydrodynamic effects 
from the effect of membrane shape deformations. 

In DPD, the interaction between any two particles, 
within a range ro, is linearly repulsive. The pairwise inter- 
action leads to full momentum conservation, which in turn 
brings out the correct fluid hydrodynamics. The force on 




Fig. 1: Image of the model film with the bulk fluid called sol- 
vent. The yellow and red particles represent the two compo- 
nents constituting the model film, while blue ones represent 
the solvent. For clarity, only a fraction of the solvent particles 
are shown. 



a particle i is given by 

™^ = E + Fg(r,„ V,,) + F«:(r,,)] , (1) 

where Vij — Vi — Vj and Vy- = v.^ — Vj. Of the three 
types of forces acting on the particles, the conservative 
force on particle i due to j is F^ = aijUj{rij)vij where 
aij is an interaction strength and Vij — Vijjrij with 
Tij = The second type of force is the dissipative 

force FJ^- = —Vi.jbj^{rij){vi^ ■■Vij)rij, where F.y is the dissi- 
pative strength for the pair (i,j). The last is the random 
force Ffj = aij{At)^^^^uj{rij)(ijrij , where cjij is the am- 
plitude of the random noise for the pair and 
is a random variable with zero mean and unit variance 
which is uncorrelated for different pairs of particles and 
different time steps. The dissipative and random forces 
act as a thermostat, provided the fluctuation-dissipation 
theorem afj ~ 2Tijk^T is satisfied (fee is Boltzmann con- 
stant and T is the thermostat temperature). The weight 
factor is chosen as i^{Tij) = 1 — Vij/r^ up to the cutoff 
radius rg and zero thereafter. The particle trajectories 
are obtained by solving eq. (H]) using the velocity- Verlet 
integrator. In the simulation, tq and m set the scales for 
length and mass, respectively, while k^T sets the energy 

1/2 

scale. The time is measured in units of r = [mrQ/k-oT^ 
The numerical value of the amplitude of the random 
force is assumed to be the same for all pairs such that 

(Tij = 3.0 [(fcBr)^m/rQ] and the fluid density is set as 
p = 3.0. We set k-^T = 1 and the integration time step is 
chosen to be = O.Olr. 

The thin film is constructed by placing particles in 
the xy-plane in the middle of the simulation box (see 
fig. [Ij. Owing to the structureless representation of the 
constituent particles, we apply an external potential so as 
to maintain the film integrity. This is done by fixing the 
z-coordinates of all the film particles. It is known that 
confinement of simulation particles between walls leads to 
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Fig. 2: The snapshots for a 70:30 mixture (yellow and red) 
undergoing phase separation at f = 0, 150 and 1000 (top to 
bottom) for a pure 2D (left column) and quasi-2D system with 
— 40 (right column). The above snapshots are from one of 
the ten independent trials that were conducted. 



a reduction in the solvent temperature near the wall |20j . 
However, since we allow for the in-plane motion of the 
film particles, the solvent temperature is found to be only 
2% less than the bulk near the 2D film. Hence we con- 
sider that this effect is negligible. Our work involves the 
systematic variation of the height of the simulation box 
starting from the pure 2D case. In the absence of solvent, 
we work with a 2D-box of dimensions Ly = 80 x 80 
with 19, 200 particles constituting the film. For the quasi- 
2D studies, we add solvent particles S above and below 
the model film, and increase the height of the box as 
Lz = 5,20 and 40. For all the cases there are 19,200 
film particles. The largest box size (L^ = 40) has 748, 800 
solvent particles. The box with height — 40 is found 
to be sufficiently large enough to prevent the finite size 
effect which affects the solvent-film interaction. The sys- 
tem is then subject to periodic boundary conditions in all 
the three directions. For phase separation simulations, we 




Fig. 3: The average domain size i? as a function of time t for a 
70 : 30 off-critical mixture. The upper black curve is the pure 
2D case showing an a = 1/2 scaling, and the lower red curve 
is the quasi- 2D case when Lz — 40 showing a distinct a = 1/3 
scaling. The inset shows the zoomed in portion with different 
box heights in the z-direction, i.e., = 0, 5, 20 and 40 starting 
from the top black curve. 



introduce two species of film particles A and B. The in- 
teraction parameter between various particles are given by 
QAA = clbb = ass = aAs = ass = 25 and gab = 50. In 
order to do a quench, the film is first equilibrated with a 
single component, following which a fraction of the parti- 
cles are instantaneously changed to the B type. 

Phase separation. — First we describe the results of 
the phase separation dynamics. The snapshots for A : B 
composition set to 70 : 30 (off-critical mixture) is shown 
in fig. [5] for both pure 2D case (left column) and quasi- 
2D case with Lz = 40 (right column). Qualitatively, it is 
seen that the domains for the quasi-2D case are smaller 
in size when compared at the same time step. We also 
monitor the average domain size R{t) which can be ob- 
tained from the total interface length L{t) between the 
two components. This is because R{t) and L(t) are re- 
lated by L(t) = 2TTN{t)R{t), where N{t) is the number of 
domains. The area occupied by the i?-component is given 
by ^ = 'nN{t)R^{t) which is a conserved quantity. Then 
we have 

R{t) = 2AlL[t). (2) 

When the domain size grows as i? ~ t", one has L ^ 
i-" and N - t'^". The domain size R{t) for 70 : 30 
mixture is shown in fig. [3l In this plot, average over 10 
independent trials has been taken. It can be seen that the 
pure 2D case has a growth exponent a = 1/2. Upon the 
addition of solvent, we observe that the exponent shifts to 
a lower value of a = 1/3. This exponent is reminiscent of 
the phase separation dynamics of an off-critical mixture in 
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a= 1/3 
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Fig. 4: The average domain size 7? as a function of time t for 
a 50 : 50 critical mixture. The upper black curve is the pure 
2D case showing an a = 1/2 scahng, and the lower red curve 
is the quasi-2D case when Lz = 40 showing a distinct a = 1/3 
scaling. The inset shows the zoomed in portion with different 
box heights in the z-direction as in fig. [S] 



3D. Upon systematically increasing the amount of solvent 
in the system by changing the height L^, we can see a 
clear deviation from the pure 2D behavior (see the inset 
of fig. 121). There is no further change if Lz is increased 
beyond 40. We note that the scaling regime covers about 
one decade in time, which is similar to that previously 
shown in the literature [T^ . A larger system size L^xLy = 
200 X 200 also produced the same scaling for the pure 2D 
case, which demonstrates that finite-size effects are small. 
However, we present here only the results for the 80 x 80 
system in 2D, because this is the system size studied for 
the quasi-2D case with a bulk fluid, which requires a large 
amount of particles in 3D. 

In fig. m we show the result for a component ratio of 
50 : 50 (critical mixture). In this case, the growth expo- 
nent for the pure 2D case is less obvious owing to rapid 
coarsening of the domains. However, by simulating a big- 
ger system 200 x 200 with the same areal density, an 
q: = 1/2 exponent is indeed obtained. Similar to the off- 
critical case, the growth of the domains is slowed down 
by the addition of solvent and the exponent is reduced to 
a = 1/3. These results indicate that solvent is responsible 
for slowing down the growth dynamics. 

The observed exponent a = 1/2 in pure 2D systems can 
be explained in terms of the Brownian coagulation mech- 
anism [l^. From dimensional analysis, the 2D diffusion 
coefficient of the domain is given by D2 k^T where 
77 is the film 2D viscosity. Using the relation 

i?2 _ _ ikBT/Tj)t, (3) 
we find R ^ t^^^. For 3D systems, on the other hand, the 



diffusion coefficient of the droplet is inversely proportional 
to its size, D3 ^ ^/R, a well-known Stokes-Einstein rela- 
tion. Hence the Brownian coagulation mechanism in 3D 
gives rise to an exponent a = 1/3. (In general, the ex- 
ponent is a = 1/d, where d is the space dimension.) The 
change in the exponent from a = 1/2 to 1/3 due to the 
addition of solvent implies the crossover from 2D to 3D 
behaviors of the phase separation dynamics even though 
the lateral coarsening takes place only within the 2D ge- 
ometry |21| . Therefore it is necessary to examine the size 
dependence of the domain diffusion coefficient in quasi-2D 
systems. This can be calculated by tracking the mean- 
squared displacement of domains of various radii. The 
equivalent information can be more efficiently obtained 
by calculating the two-particle longitudinal coupling dif- 
fusion coefficient in a single component film rather than 
in a binary system. This is described in the next section. 

Correlated diffusion. — Consider a pair of particles 
separated by a 2D vector s, undergoing diffusion in the 
liquid film. The two-particle mean squared displacement 
is given by [H] 



(4) 



where AsJ' is the displacement of the particle k{= 1,2) 
along the axis i{= x,y), Df^ is the diffusion tensor giving 
self-diffusion when k = I and the coupling between them 
when k I. The a;- axis is defined along the line connecting 
a pair of particles 1 and 2, i.e., s = sx. Hence, we have 
^xy = by symmetry. The longitudinal coupling diffusion 
coefficient, Df^{s) = £);J,^(si;), gives the coupled diffusion 
along the line of centers of the particles. We first describe 
the analytical expression of -D£(s) based on the Saffman 
and Delbriick (SD) theory which was originally developed 
for protein diffusion in membranes |15| . 

Since the calculation of the diffusion coefficient in a pure 
2D system is intractable due to Stokes paradox, SD cir- 
cumvented this problem by taking into account the pres- 
ence of the solvent with 3D viscosity i]s above and below 
the membrane. Suppose a point force f directed along the 
plane of the film lying in the ccy-plane acts at the origin. 
Then we seek for the velocity v(s) induced at the position 
s. According to the SD theory, it is given in Fourier space, 
v[q] = / dse~*i-^v(s), as dSHmHSI 



1 



(5) 



where G^^ is the 2D ffim analog of the Osecn tensor. In 
the above, the SD length is defined by i^^^ = v/'^Vs- 

For over-damped dynamics, we can use the Einstein re- 
lation to relate the diffusion tensor Z)£ to G^J^ [22] . After 
converting Gfj'^[q] into real space, we obtain 



k^TGl^isx) 



keT 



7rHi(i^s) ttYi(vs) 



vs 



{vsf 



(6) 
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Fig. 5: Longitudinal coupling diffusion as a function of 
particle separation s. The upper circles are data for the pure 
2D case. The lower squares correspond to the case with solvent 
when Lz — 40. The upper solid line is the fit by eq. 0, and 
the lower solid line is the fit by eq. ((6|. The dashed line shows 
the 1/s dependence. 



where Hi and Yi are Struve function and Bessel function 
of the second kind, respectively. At short distances s ^ 
, the asymptotic form of the above expression becomes 



Dlis) 



kBT 



(7) 



where 7 = 0.5772 • • • is Euler's constant. At large inter- 
particle separations s ^ on the other hand, eq. © 
reduces to 



(8) 



2'Krjvs iTTrjsS 

showing the asymptotic 1/s decay which reflects the 3D 
nature of this limit. Notice that eq. (|5]) depends only on 
the solvent viscosity ?7s but not on the film viscosity 77 any 
more. 

In fig. [5l we plot the measured longitudinal coupling 
diffusion coefficient Df^ as a function of 2D distance s. In 
these simulations, we have worked with only single compo- 
nent films with the same system sizes and number of par- 
ticles as those used for the phase separation simulations. 
We have also taken average over 20 independent trials. In 
the pure 2D case without any solvent, shows a loga- 
rithmic dependence on s. This is consistent with eq. (0 
obtained when the coupling between the film and solvent 
is very weak so that the film can be regarded almost as a 
pure 2D system. Using eq. (O as an approximate expres- 
sion, we get from the fitting as kBT/Amj w 0.89 x 10^^ 
and ~ 20. In an ideal case, the SD length should 
diverge due to the absence of solvent. The obtained finite 
value for i^'^ is roughly set by the half of the system size 
in the simulation. 



When we add solvent (L^ = 40), the £)£ is decreased 
and no longer behaves logarithmically. In this case, we 
use the full expression eq. (O for the fitting, and obtained 
kBT/Airrj « 1.35 x 10~^ and v^^ w 1. In the above two 
fits we have neglected the first two points as they lie out- 
side the range of validity, s :§> 1, of eq. ([H]) [Hj. Since 
v"^ « 1 when the solvent is present, the data shown in 
fig. [S] are in the crossover region, s > v~^, showing an 
approach towards the asymptotic 1/s dependence as in 
eq. ([8]). Hence we conclude that the solvent brings in the 
3D hydrodynamic property to the diffusion in films. This 
is the reason for the 3D exponent a = 1/3 in the phase 
separation dynamics, and justifies that it is mainly driven 
by Brownian coagulation mechanism. 

In our simulations the film and the solvent have very 
similar viscosities. This sets the SD length scale to be 
of the order of unity, which is consistent with the value 
i^^^ « 1 obtained from the fitting. As explained above, 
the fit also provides the 2D film viscosity as 77 « 6, and 
hence we obtain as 77s ~ 3. This value is in reasonable 
agreement with the value r/g « 1 calculated in ref. [Mj by 
using the reverse Poiseuillc flow method. The reason for 
the slightly higher value of ?7s in our simulations is that the 
tracer particles are of the same size as the film particles. 
This may lead to an underestimation of the correlated dif- 
fusion coefficient. In real biomcmbranes sandwiched by 
water, however, the value of the SD length is much larger 
than the lipid size, and is in the order of micron scale |15| . 
Hence the 3D nature of hydrodynamics can be observed 
only for large enough domains observed under optical mi- 
croscopes [5]. 

Discussion. — Several points merit further discus- 
sion. The growth exponents obtained from our simula- 
tions for critical mixtures are the same for the off-critical 
case, namely a = 1/2 without solvent and a = 1/3 with 
solvent. A previous DPD study by Laradji and Kumar on 
phase separation dynamics of two-component membranes 
(both critical and off-critical cases) used a coarse-grained 
model for the membrane lipids |18j . In their model, the 
self-assembly of the bilayer in the presence of solvent is 
naturally taken into account. The exponent for the off- 
critical case a = 1/3 is the same as that obtained in our 
study, although they attributed this value to the LS mech- 
anism. For critical mixtures in the presence of solvent, 
they obtained a different value a = 1/2. A suitable expla- 
nation for this exponent was not given in their paper. 

As a related experimental work, the diffusion of tracer 
particles embedded in a soap film was recently re- 
ported When the diameter of the tracer particles 
is close to the thickness of the soap film, the system shows 
a 2D behavior. On the other hand, if the particle diameter 
is much smaller than the soap film thickness, it executes 
a 3D motion. On systematically increasing the soap film 
thickness, they identified a transition from 2D to 3D na- 
ture. In this paper, we have demonstrated the analogue 
for a 2D liquid film-solvent system using DPD simulations. 
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In the SD theory, the bulk fluid is assumed to occupy 
an infinite space above and below the membrane. The 
situation is altered when the solvent and the membrane 
are confined between two solid walls [22- If the distance 
H between the membrane and the wall is small enough, 
we arc in a situation similar to that described by Evans 
and Sackmann [55]. The Oseen tensor G^^ in this case is 
defined through the relation [37] , 



1 



5 \ f- 



(9) 



where the new length scale is defined as = y''j]H/2j]s- 
Notice that is the geometric mean of and H [^ . 
Following the same procedure as in the previous section, 
the longitudinal coupling diffusion coefficient can be ob- 
tained as 
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where Ki is modified Bessel function of the second kind. 
At short distances s <C k^^, we have 



ln( — 

KS 



(11) 



which is almost identical to eq. (O except v is replaced 
now by k. At long distances s S> on the other hand, 
we get 

kBT k^TH 



Dlis) 



27r?7K2 



(12) 



which exhibits a l/s^ dependence. This is in contrast to 
eq. ([8]). Following the similar scaling argument as before, 
we predict that, in the presence of solid walls, the domain 
growth exponent should be a = 1/4 within the Brownian 
coagulation mechanism. In biological systems, the above 
model with solid walls can be relevant because the cell 
membranes are strongly anchored to the underlying cy- 
toskeleton, or are tightly adhered to other cells. 

In summary, we have shown that the bulk fluid has a 
significant effect on the phase separation dynamics of thin 
liquid films through a simple quasi-2D simulation model. 
We have demonstrated the change in the growth expo- 
nents from 2D to 3D by the addition of bulk fluid. This 
is further justifled by the two-particle correlation studies, 
which showed that the longitudinal coupling diffusion co- 
efficient changes from a logarithmic dependence for the 
pure 2D case to an algebraic one for the quasi-2D case. 
Future directions include investigating the role of out-of- 
plane fluctuations and the effect of boundary walls on the 
phase separation. 
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